Segregation process and phase transition in cyclic predator-prey models 

with even number of species 
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We study a spatial cyclic predator-prey model with an even number of species (for n — 4, 6, and 8) that 
allows the formation of two defective alliances consisting of the even and odd label species. The species are 
distributed on the sites of a square lattice. The evolution of spatial distribution is governed by iteration of two 
elementary processes on neighboring sites chosen randomly: if the sites are occupied by a predator-prey pair 
then the predator invades the prey's site; otherwise the species exchange their site with a probability X. For low 
X values a self-organizing pattern is maintained by cyclic invasions. If X exceeds a threshold value then two 
types of domains grow up that formed by the odd and even label species, respectively. Monte Carlo simulations 
indicate the blocking of this segregation process within a range of X forn = 8. 

PACS numbers: 87.23.Cc, 89.75.Fb, 05.50.-l-q 



Interesting phenomena in the cyclic predator-prey systems 
initiated a progressive research in the last decades. In the sim- 
plest case the system has three species dominating cyclically 
each other, that is. Si beats 5*2 beats S3 beats 5*1. Cyclic dom- 
inance occurs in many systems including biological models 
iilSlllH] and evolutionary games |5, 6]. It is well known that 
such a cyclic dominance can sustain all of the three species in 
a well-mixed community iHIl]- In the spatial version of this 
system the individuals stay on the sites of a lattice [9] and 
the cyclic invasions maintain a self-organizing pattern in the 
spatial distribution of species. Nowadays, this system is fre- 
quently referred to as an evolutionary spatial Rock-Scissors- 
Paper game (for a review see [10.1 ). 

Within the framework of mean-field approximation (assum- 
ing well-mixed population) the species densities can be either 
stationary or oscillatory (for the latter case both the sum and 
product of species densities are conserved quantities) Isi lTlll . 
Furthermore, these systems exhibit unusual response to the 
variation of invasion rates or to any external support pl^ fTsll . 
For finite population the system evolves into one of the ho- 
mogeneous state via a Moran process ITii ITsll . On the other 
hand, when the species are distributed on a lattice the cyclic 
invasions yield a self-organizing pattern lfTl [T6ll providing sta- 
bility against some types of external invaders jlTi [Tsl] if the 
spatial dimension is larger than one (d > 1). For the one- 
dimensional lattice the dynamical rule results in growing do- 
mains and super-domains |9, 19, 20]. On the Bethe lattice one 
can observe limit cycle behavior (for a degree of z = 3 or 
4) and the oscillation grows until the system reaches one of 
the homogeneous states if z > 6 ll2Tll22ll . More complicated 
behavior is reported on some types of directed graphs |23]. 

Starting from the simplest model there are many ways for 
the generalization of the evolutionary Rock-Scissors-Paper 
game. A straightforward possibility is to increase the number 
n of species (states) so that the cyclic invasion remains valid, 
i.e.. Si invades S2 invades S3 etc. and finally Sn invades 
Si . In the one-dimensional lattice domain growth can be ob- 
served until n < 5, otherwise the spatial distribution tends to- 
wards a frozen domain structure where the species staying in 
neighboring domains cannot invade each other (shortly, they 



are neutral). Similar fixation process was reported for both 
the two- and three-dimensional lattices if n exceeds a criti- 
cal value dependent on d Further relevant observation 
is related to the parity of n affecting some features (e.g., the 
sensitivity to the inhomogeneous invasion rates) in the system 
ifTTl fisll . Consequences of supplementary microscopic pro- 
cesses are also investigated. For example, one can introduce 
mutation, extinction, and empty sites to which the neighboring 
species can jump, allowing spatial mixing ll26ll27ll . Now our 
efforts will be concentrated on a lattice Lotka-VolteiTa sys- 
tem with an even number n of species dominating cyclically 
each other while the stochastic local mixing is described by 
a site exchange mechanism between neutral pairs. The sim- 
plest model where such a process can be introduced is the 
four-species version that exhibits a phase transition when in- 
creasing the strength of mixing |26, 28]. More precisely, for 
low rate of mixing one can observe a self-organizing pattern 
resembling the sample of evolutionary Rock-Scissors-Paper 
game. On the contrary, if the strength of mixing exceeds a 
threshold value then phase segregation occurs, i.e., the well- 
mixed odd end even label species form growing domains and 
finally one of these states dominates the whole finite system. 
Our primary interest was to investigate how the critical value 
of mixing decreases when increasing the number of species. 
During this study, however, an unexpected intermediate phase 
was found for n = 8 and 10 as will be reported below. 

We consider a lattice Lotka-VolteiTa model on a square 
lattice where each site x is occupied by a single individual 
belonging to one of the n species (n is even), that is, the 
distribution of species can be described by the site variables 
Sx — 1, . . . ,n referring to the label of species. The time evo- 
lution of the species distribution is determined by invasions 
between the site x (chosen at random) and one of the (ran- 
domly chosen) neighboring sites y if the sites are occupied by 
a predator-prey pair, i.e., the {s^, Sy) [and also the {sy, s^)] 
pair transforms into {sx, s^) if the species is the predator 
of Sy. The predator-prey relation is defined by a cyclic food 
web and the invasion rates between any predator-prey pair are 
equivalent and chosen to be unity. Besides it, if the species 
Sx and Sy are neutral then the species may exchange their site 
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with a probability X characterizing the strength of mixing. 
Obviously, nothing happens if ~ Sy. During the time unit 
(called Monte Carlo step, shortly MCS) the above elementary 
process is repeated once on average for each site. 

For the Monte Carlo (MC) simulations the system is started 
from a random initial state (providing equivalent (average) 
number of individuals for the species) on a square lattice (con- 
sisting of = Ly.L sites) with periodic boundary conditions. 
When repeating the above defined elementary step (invasion 
or local mixing) for sufficiently large size the system evolves 
into a stationary state that can be described by the average 
density pi of species i (satisfying the condition X]r=i Pi ~ 
and by the pair configuration probabilities p2 {hj) of finding 
species i and j on two neighboring sites. In the limit N ^ oo 
each species density remains constant, that is pi{t) — 1/n. 
Whereas the ordering processes can be well quantified by con- 
sidering the time-dependence of pair configuration probabil- 
ities. Due to the symmetries we will distinguish two basic 
quantities. The predator-prey pak probability is given as 



Ppp 



(1) 



where the time-dependence of p2 is not denoted and i + 1 
means 1 for i = n. The other analogous quantity is the neutral 
pair probability that can be expressed as 



(2) 



i=l fe=2 



where z + A; is cyclically reduced to the range 1 to n. For the 
MC simulations these quantities are determined by averaging 
over a suitable sampling interval. 

For small size (e.g., L < 10) the system quickly develops 
into one of the states in which several species and simultane- 
ously the predator-prey invasions are missing, thus the compo- 
sition remains constant. The average fixation (transient) time 
increases fast with the linear size L. For sufficiently large size, 
however, most of the latter phases can also be observed locally 
within small patches and the evolution of spatial distribution 
is affected by the competition between these phases lfloll28|] . 

First we briefly recall the MC results obtained for n = 
4 ll28ll . According to the simulations all the four species 
are sustained by the cyclic invasions if X < ^c(4) = 
0.02662(2). Within this region of X the value of p„ in- 
creases monotonously with X meanwhile an opposite trend 
occurs in ppp as demonstrated in Fig. [T] At the critical value 
of X a sudden change occurs in both quantities because, 
through a domain growing process, the finite system evolves 
into a state where either the odd or the even labelled species 
form a well-mixed phase, that is, pi(oo) ~ Psioo) — 1/2 
and p2(oo) ~ P4{oo) = or pi(oo) = P3(oo) = and 
P2(oo) = p4(oo) — 1/2. In both phases (cx)) = and 
p„(oo) = l/2(ifX>X,(4)). 

To support our previous MC results we have also performed 
generalized mean-field approximations at different levels. In 
the case of the four-species version of this model, all the con- 
figuration probabilities on 2 x 1-, 2 x 2-, and 3 x 3-site clus- 
ters have been evaluated (for details of the method see ifToll ). 
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FIG. 1 : MC data for the predator-prey pair ppp (diamonds) and neu- 
tral pair pn (squares) probabilities as a function of X in the final 
stationary state for n = 4. Dotted, dashed, and solid lines represent 
the prediction of generalized mean-field approximation for ppp at the 
levels of 2 x 1-, 2 x 2-, and 3 x 3-site clusters. The left arrow in- 
dicates the position of Xc, the right arrow shows the corresponding 
prediction obtained by the 3 x 3-site approximation. 



Evidently, the traditional mean-field approach (one-site ap- 
proximation) cannot take into account the effect of mixing. 
At the levels of 2 x 1- and 2 x 2-site cluster approximation 
this method is not capable to describe the transition observed 
by MC simulations although the solution of many well-mixed 
phases exists. The more accurate 3 x 3-site approximation pre- 
dicts a phase transition at X = X^^^'(4) = 0.1285(5). The 
large deviation from the MC results indicates the importance 
of consecutive elementary steps yielding important correla- 
tions in the spatial distribution. 

For sufficiently strong mixing the phase segregation pro- 
cess seems to be a robust phenomenon as it is observed for 
other dynamical rules lfToll26ll27ll . It turned out that the well- 
mixed phases of the odd (as well as even) label species can be 
considered as a defensive alliance because within this spatial 
association the species guard each other against the external 
invaders. For example, if species 1 is attacked by an individ- 
ual of species 4 then one of the neighboring species 3 strikes 
back within a short time. Due to the cyclic symmetry species 
1 guards species 3 against species 2 and similar mechanism 
protects the well mixed spatial association of the even label 
species. 

One can easily check that the above concept of defen- 
sive alliances remains valid for any even number of species. 
Namely, within the well-mixed phase of the odd label species 
the species protect each other against the invasion of even 
label species and vice versa. Thus for larger n one can ex- 
pect a similar phase transition from the cyclic self-organizing 
pattern to the phase segregation phenomenon if we increase 
the value of X. MC simulations confirm this expectation for 
n = 6 as illustrated in Fig. |2] The predator-prey probability 
(Ppp) drops suddenly to zero at X = Xc{6) = 0.00654(2). 
The generalized mean-field analysis of this system is not per- 




FIG. 2: The non-vanishing predator-prey pair probabilities vs. X in 
the stationary states forn = 4 (triangles), 6 (squares), and 8 (circles). 

Notice that many possible solutions emerge if n is in- 
creased. For example, the well-mixed phases of the mutually 
neutral species (e.g., 1+3+6 for n = 8) with arbitrary compo- 
sition are stationary states. All these states can occur in the 
MC simulations for small sizes (L < 10) and can also be re- 
produced by the generalized mean-field methods. Indeed, the 
solutions of the subsystems (several species are missing) are 
solutions for the whole system too. Despite of the large num- 
ber of possible solutions, the visualization of spatial distribu- 
tions in the MC simulation indicates the presence of only the 
above mentioned relevant phases for sufficiently long times. 

Surprisingly, for the eight-species system three types of 
phases (behaviors) can be distinguished as illustrated in Fig.|2] 
These MC data are obtained for L — 400 or 600. For these 
sizes the domain growing process ends within a few million 
MCS and one of the four-species defensive alliances (consist- 
ing of only the odd or even label species) prevails the whole 
system in the final state if X > Xc2{8) ~ 0.0042(5). The 
self-organizing spatio-temporal pattern can be maintained by 
the cycHc invasions until X < Xci{S) ~ 0.00285(3). As a 
consequence, the MC simulations reveal the appearance of an 
intermediate phase within a region of Xci(8) < X < Xc2{8) 
where the domain growing process stops (or becomes ex- 
tremely slow). 

To illustrate the formation of defensive alliance in the in- 
termediate region (Xd < X < Xc2) the odd and even label 
species are denoted by light and dark colors in Fig. [3] The 
snapshot for t = 400 000 MCS shows clearly that the territo- 
ries of defensive alliances are separated by a boundary layer 
where the cyclic invasions govern the time evolution. Similar 
patterns (with a thickness dependent on X) can be observed 
during the domain growing process for X > Xc2{8). We 
have to emphasize that these boundary layers play crucial role 



FIG. 3: (Color online) Spatial distribution of species after 400 000 
MCS for X = 0.003 if initially the eight species were distributed 
randomly on a square lattice. The cyclic dominance between the 
eight (colored) species is indicated at the top and the snapshot shows 
a 400 X 400 portion of the whole system with a size of 1600 x 1600. 



in the formation of final pattern. Namely, these layers serve as 
a symmetric species reservoir for both defensive alliances and 
help the equalization of their composition via diffusion. 

As mentioned above the domain growing process can be 
investigated quantitatively by recording the probability of 
predator-prey pairs because such a constellation occurs ex- 
clusively within the boundary layers. One can think that the 
inverse of Ppp is proportional to the average linear size of do- 
mains of defensive alliances. 

Figure |4] shows some typical behaviors when considering 
the time-dependence of ppp for such a system size where 
L is significantly larger than the average domain size at the 
end of simulation (here at t = 10^ MCS). In order to sup- 
press the short-time fluctuations, the numerical data of Ppp{t) 
are averaged over a time interval with a typical width of 
t^, ~ min(i/20, 1000 MCS). The upper curve (X = 0.0025) 
in Fig. m illustrates that the frequency of invasions reaches a 
high stationary value (ppp ~ 0.11) characterizing the cyclic 
self-organizing pattern if X < Xci{8). The lowest curve 
(for X — 0.008) indicates a typical domain growing process 
when the asymptotic behavior (ppp ^ t^^^'^) becomes similar 
to phase ordering process with non-conserved dynamics ll29tl . 
From the plotted MC data at X = 0.005 and 0.008 one can 
suggest similar asymptotic behavior Data for X = 0.003, 
however, indicates that the domain growth is stopped and the 
process itself is resembling the segregation of a water-oil mix- 
ture in the presence of surfactant |30, 31]. Notice that in this 
spatio-temporal pattern all the eight species remain alive, thus 
the spontaneous formation of this inhomogeneous pattern ex- 



4 




t [MCS] 



FIG. 4: Log-log plot for the time-dependence of the predator-prey 
pair probabilities in the eight-species system for X — 0.0025, 0.003, 
0.004, 0.0045, 0.005, and 0.008 (from top to bottom). The MC re- 
sults are obtained on a square lattice with a linear size of L = 2800 
and the plotted data are smoothed by averaging over a time-window 
(typically At ~ t/40 maximum 1000 MCS). Dashed line shows the 
slope of —1/2 characterizing domain growth driven by the decrease 
of interfacial energy 1 29ll . 

emplifies a way how the biodiversity can be maintained for 
a long time. We have to emphasize that, in the absence of a 
clear theoretical explanation of this phenomenon, we cannot 
exclude the appearance of a slower domain growing process 
for a longer time scale {t > 10^ MCS). 

Within the intermediate region of X the reproduction of 
numerical data is poor (see Fig. |2]i because the system behav- 
ior is perturbed by large and slow fluctuations (see data for 
X — 0.003 in Fig.|4]i. To have a deeper insight into this order- 
ing process the function ppp{t) is plotted in Fig.|5]for different 
system sizes. 

Figures |4] and |5] indicate clearly that the sufficiently large 
domains of the well mixed phases of odd and even label 
species are formed after a transient time as long as 20,000 
MCS. The symmetric composition of these domains (namely, 
Po — P2 — Pi = P6 = 1 /4) is ensured by the long contact 
(interaction) with the boundary layers serving as symmetric 
species reservoirs |26]. At the same time, the final spatial 
structure of the boundary layers is also affected by their inter- 
actions with the well-mixed phases of neutral species. These 
mutual effects can be a cause of the blocking of domain grow- 
ing process in the intermediate region of X. The mentioned 
process is observed if the system size L exceeds significantly 
the typical domain size (/ ~ 100 lattice unit in Fig.O. In the 
opposite case this proper structure cannot build up because 
the system develops into a state consisting of two, three or 
four neutral species. For example, the left curve {L = 10) in 
Fig.|5]represents a fast evolution into a final state where about 
one third (two percent) of runs end with three (two) neutral 
species. Although the probability of finding four neutral sur- 
viving species increases with system size, the composition of 



the final state is far from being symmetric if I ^ L. 
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FIG. 5: MC data of Ppp(t) for X = 0.003 when varying the system 
size. Thick lines from left to right show the results for L — 10, 30, 
and 100 after averaging over 10", 10^ and 10^ runs, respectively. To 
demonstrate the large fluctuations in the fixation time the subsequent 
three dashed lines illustrate the results of three runs for L — 200. The 
last thick line (representing a single run for L = 1000) resembles 
data (obtained for L = 2800) plotted in Fig.|4] 



Due to the extremely long transient times at the boundaries 
of the intermediate region, more accurate determination of the 
critical values and systematic analysis of the corresponding 
phase transitions exceed our computing capacity. 

The preliminary results indicate similar behavior for n — 
10. In this case the system exhibits longer relaxation (re- 
lated to the slower formation of the corresponding defensive 
alliances) making the rigorous analysis more difficult. 

In summary, the present work is focused on the spatial for- 
mation of two defensive alliances on a two-dimensional lattice 
Lotka- Volterra model with even number n of species invading 
cyclically each other with the same rate. The introduction of 
local mixing (with strength characterized by the site exchange 
probability X between the neutral species residing on neigh- 
boring sites) supports the formation of the well-mixed distri- 
bution of odd or even label species representing two equiv- 
alent defensive alliances. Phase segregation process is ob- 
served if the mixing rate exceeds a threshold value dependent 
on n. According to our MC simulations this system evolves 
into a pattern (for n > 6) where the domains of defensive 
alliances are separated by boundary layers having a different 
structure. This phenomenon raises further questions about the 
role of boundary layers in these types of complex systems. 



Acknowledgments 

This work was supported by the Hungarian National Re- 
search Fund (Grant No. T-47003). 



5 



[1] B. Kerr, M. A. Riley, M. W. Feldman, and B. J. M. Bohannan, 

Nature 418, 171 (2002). 
[2] R. Durrett and S. Levin, Theor. Pop. Biol. 53, 30 (1998). 
[3] C. R. Johnson and I. Seinen, Proc. Roy. Soc. Lond. B 269, 655 

(2002). 

[4] M. Mobilia, I. T. Georgiev, and U. C. Tauber, J. Stat. Phys. 128, 
447 (2007). 

[5] M. Nowak, J. Theor. Biol. 142, 237 (1990). 

[6] A. Traulsen and J. C. Claussen, Phys. Rev. E 70, 046128 (2004). 

[7] R. M. May and W. J. Leonard, SIAM J. Appl. Math. 29, 243 

(1975). 

[8] J. Hofbauer and K. Sigmund, Evolutionary Games and Pop- 
ulation Dynamics (Cambridge University Press, Cambridge, 
1998). 

[9] K. Tainaka, Phys. Rev. Lett. 63, 2688 (1989). 
[10] G. Szabo and G. Path, Phys. Rep. in press (2007). 
[11] Y. Itoh, Progr. Theor. Phys. 78, 507 (1987). 
[12] K. Tainaka, Phys. Lett. A 176, 303 (1993). 
[13] M. Frean and E. D. Abraham, Proc. R. Soc. Lond. B 268, 1 
(2001). 

[14] Y. Itoh and K. Tainaka, Phys. Lett. A 189, 37 (1994). 

[15] M. Mobiha, I. T. Georgiev, and U. C. Tauber, Phys. Rev. E 73, 

040903(R) (2006). 
[16] K. Tainaka, Phys. Rev. E 50, 3401 (1994). 



[17] M. C. Boerlijst and P Hogeweg, Physica D 48, 17 (1991). 
[18] G. Szabo and T. Czaran, Phys. Rev. E 63, 061904 (2001). 
[19] L. Frachebourg, P. L. Krapivsky, and E. Ben-Naim, Phys. Rev. 

Lett. 77, 2125 (1996). 
[20] L. Frachebourg, P. L. Krapivsky, and E. Ben-Naim, Phys. Rev. 

E 54, 6186 (1996). 
[21] K. Sato, N. Konno, and T. Yamaguchi, Mem. Muroran Inst. 

Tech. 47, 109 (1997). 
[22] A. Szolnoki and G. Szab6, Phys. Rev. E 70, 027101 (2004). 
[23] C.-Y. Ying, D.-Y. Hua, and L.-YWang, J. Phys. A: Math. Theor. 

40, 4477 (2007). 

[24] L. Frachebourg and P L. Krapivsky, J. Phys. A 31, L287 (1998). 
[25] K. Sato, N. Yoshida, and N. Konno, Appl. Math. Comp. 126, 
255 (2002). 

[26] G. Szabo and G. A. Sznaider, Phys. Rev. E 69, 03I9I1 (2004). 
[27] M. He, Y. Cai, Z. Wang, and Q.-H. Pan, Int. J. Mod. Phys. C 

16, 1861 (2005). 
[28] G. Szabo, J. Phys. A: Math. Gen. 38, 6689 (2005). 
[29] A. J. Bray, Adv. Phys. 43, 357 (1994). 

[30] G. Gompper and M. Schick, Self-assembling amphiphilic sys- 
tems (Academic Press, London, 1994). 

[31] J. R. Henriksen, M. C. Sabra, and O. G. Mouritsen, Phys. Rev. 
E 62, 7070 (2000). 



